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ABSTRACT 

We revisit our original papers on the burst mode of accretion by incorporating a detailed energy 
balance equation into a thin-disk model for the formation and evolution of circumstellar disks around 
low-mass protostars. Our model includes the effect of radiative cooling, viscous and shock heating, 
and heating due to stellar and background irradiation. Following the collapse from the prestellar 
phase allows us to model the early embedded phase of disk formation and evolution. During this 
time, the disk is susceptible to fragmentation, depending upon the properties of the initial prestellar 
core. Globally, we find that higher initial core angular momentum and mass content favors more 
fragmentation, but higher levels of background radiation can moderate the tendency to fragment. 
A higher rate of mass infall onto the disk than that onto the star is a necessary but not sufficient 
condition for disk fragmentation. More locally, both the Toomre Q-parameter needs to be below a 
critical value and the local cooling time needs to be shorter than a few times the local dynamical 
time. Fragments that form during the early embedded phase tend to be driven into the inner disk 
regions, and likely trigger mass accretion and luminosity bursts that are similar in magnitude to FU- 
Orionis-type or EX-Lupi-like events. Disk accretion is shown to be an intrinsically variable process, 
thanks to disk fragmentation, nonaxisymmetric structure, and the effect of gravitational torques. The 
additional effect of a generic a- type viscosity acts to reduce burst frequency and accretion variability, 
and is likely to not be viable for values of a significantly greater than 0.01. 

Subject headings: accretion, accretion disks — hydrodynamics — instabilities — ISM: clouds — stars: for- 
mation 



1. INTRODUCTION 

Typical rotation rates of ~ 1 km s~^ p c~^ ^ 10~^^ rad 
s~^ m easured in molecula r cloud cores (|Goodman et al.l 
[1993; Caselli et al. 2002) are sufficient to provide a 
significant angular momentum barrier to star formation. 
Most of the infalling matter will land on a protostellar 
disk rather than directly onto a protostar, since mag- 
netic braking is rendered ineffective by ohmic dissipation 
in the near-stellar environment. Therefore, the early 
phase of disk formation and evolution holds the key 
to understanding stellar mass accumulation, and sets 
the initial conditions for a later stage of disk evolution 
during which planets may form by core accretion 
(|Lissauerl Il993f ). The early disk phase is character- 
ized by episodic accret ion, as predicted theor e tically 
in our earlier papers (|Vorobvov fc Basul I2005L l2006f ) 
and inferred observationally by compilin g; luminosity 
distributions of young stellar objects (jEnoch et al.l 
[2009). Furthermore, the FU Orionis stars, named 
after the prototype FU Ori (Herbig 1977), provide 
direct evidence of transient luminosity variations (3 — 6 
mag during < 100 yr). These luminosity bursts have 
been associated with a sharp i ncrease of the mass 
accretion rate onto the protostar (Hartma nn fc KenyonI 
|l996) and various physical mechanisms have been 

proposed to ex pl ain t his phe nomenon (see e.g., 
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"Beh fc Lin>'1994«; "Armitage et al.' 2001; 'Lodato fc Clarkd 
2004: Vorobvoy fc Basu 2005, 2006; Pfalzner et al. 20W 
IZhu et al.ll2009l . 120101 : iForgan fc Rice 2010; Bate 20l3k 
Ou r earlier calculations ( Vorobyov fc Basuii2005u i2QQ6l 
l2007l . [2008) have revealed the importance of studying 
disk evolution using a self-consistent method of follow- 
ing the collapse of an initial prestellar core. Disk for- 
mation occurs after a central stellar core has formed, 
but the disk continues to gain mass from the surround- 
ing infalling envelope. Under certain conditions, this 
leads to disk fragmentation and the development of the 
"burst mode" of accretion, during which fragments are 
driven onto the protostar and episodic high mass accre- 
tion events actually account for the majority of mass ac- 
cumulation onto the pro tosta r. For example, F ig. 1 of 
IVorobv ov fc Basul (|2007 ^ and 'Vorobvovl (|2009bf ) reveals 
the correlation of the burst mode with infall from the 
envelope onto the disk. The above models were charac- 
terized by a large dynamic range of spatial and tempo- 
ral scales, so that core collapse from ~ 10^ AU scales 
down to an inner sink cell of size 5 — 10 AU was re- 
solved, and the evolution followed for up to several Myr 
after the formation of a central protostar. Aside from 
the initial discovery of the burst mode, long-term evo- 
lution revealed that, even after the burst mode ceases, 
the disk settles into a self-regulated mode in which the 
Toomre-Q parameter stays near the critical value. In 
this phase, residual accretion due to gravitational torques 
(resulting from persistent low-amplitude nonaxisymmet- 
ric structure driven by the swing amplifier effect) occurs 
at a rate that can explain observed T Tauri star mass 
accretion (Vorobyov fc Basu 2007, 2008), but with late 
time disk masses that are about an order of magnitude 
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greater than observational estimates (|Vorobvovl l2QQ9al ) 
which may re ahy represent lower hmits to the actual val- 
ues (see e.g.. [Andrews fc Williamsll2QQ7[ ). The effect of 
additional angular momen tum transport du e to an a- 
viscosity was explored by I Vorobyov fc Basul (|2QQ9al lbl). 
with a finding that such an a would likely lie in the 
range 10~^ — 10~^ to satisfy observational and theoreti- 
cal constraints. An a-viscosity in this range could begin 
to dominate gravitational torques only during the late 
evolution (Class II, or T Tauri phase) and yield a late 
time accretion rate that was a factor 2 — 3 greater than 
that due to gravitational torques alone. Values of a well 
above this range were found to lead to a lack of accretion 
variability in the early stages, and quickly lead to very 
low mass disks, but with disk lifetimes < 1 Myr that may 
be too short. 

The above calculations employed a barotropic equa- 
tion of state, which captured the basic features of 
the transition from isothermal to polytropic evolu- 
tion at high densities (at number density n > 10^^ 
cm~^) , as seen in spherical radiative transfer calcula- 
tions (|Masunaga fc Inutsukal 12000). However, such cal- 
culations could not capture the detailed thermodynam- 
ics in the vicinity of forming clumps. That can be a 
crucial effect in the development of a clump, including 
determining whether it can even form at all (G ammie 
[2QQTI : iRice at al.|[2QQ3l : iMeua et al.|[2QQ5t iCai et al.i r2QQ8) . 
Hence, a criticism of the above modeling has been that 
the clump formation and consequent burst mode may not 
be a robust result in the case of more realistic thermo- 
dynamics. We note that Vorobyov fc Basu (2006) were 
aware of this difficulty, and tested out models with high 
values of polytropic index such that the temperatures 
for densities n > 10^^ exceeded that found in radiative 
transfer calculations. The clump formation and bursts 
still occurred in those cases, thanks to significant forcing 
by mass accretion onto the disk during the early phases, 
although their number and frequency could be strongly 
reduced. Bursts were found to be robust in the context 
of those models, but their frequency depended strongly 
on thermal evolution as well as the mass and angular 
momentum content of infalling material. 

In this paper, we have made a major improvement 
by implementing the energy balance equation, which in- 
cludes radiative cooling, viscous and shock heating, and 
heating due to stellar and background irradiation. We 
continue to include angular momentum transport due to 
a generic a- viscosity term, since mechanisms other than 
gravitational torques may be at work. Numerical resolu- 
tion is also extended to greater values than in our orig- 
inal p apers on the burst mode (|Vorobvov fc Basul I2005L 
l2006f ). An important question is: will the existence of 
the burst mode be robust under these circumstances, and 
what will be its properties? We address these questions 
in Sections 3 - 5. A description of our model is in Sec- 
tion 2, and we provide extended discussion of the model 
features in Section 6. A summary of results is in Section 
7. 

2. DESCRIPTION OF THE NUMERICAL MODEL 

The main concept s of our numerical appro ach are ex- 
plained in detail in IVorob vov fc Basul (|2006f ). Here we 
review some main properties and focus mainly on the 
implementation of radiative cooling, viscous and shock 



heating, and heating due to stellar and background irra- 
diation. 

We start our numerical integration in the pre-stellar 
phase, which is characterized by a collapsing starless 
cloud core, continue into the embedded phase of star 
formation (hereafter, EPSF), during which a star, disk, 
and envelope are formed, and terminate our simulations 
in the T Tauri phase, when most of the envelope has ac- 
creted onto the forming star/disk system. In the EPSF, 
the disk occupies the innermost region of our numerical 
grid, while the larger outer part of the grid is taken up 
by the infalling envelope, which is a remnant of the par- 
ent cloud core. This ensures that the protostellar disk is 
not isolated in the EPSF but is subject to intense mass 
loading from the envelope. In addition, the mass accre- 
tion rate onto the disk Mq^v is not a free parameter of 
the model but is self-consistently determined by the gas 
dynamics in the envelope. 

We introduce a "sink cell" at r^c = 6 AU and impose a 
free infiow inner boundary condition. We monitor the gas 
surface density in the sink cell and when its value exceeds 
a critical value for the transition from isothermal to adi- 
abatic evolution, we introduce a central point-mass star. 
In the subsequent evolution, 90% of the gas that crosses 
the inner boundary is assumed to land onto the central 
star plus the inner axisymmetric disk at r < 6 AU. This 
inner disk is dynamically inactive, it contributes only to 
the total gravitational potential and is used to secure a 
smooth behavior of the gravity force down to the stellar 
surface. The other 10% of the accreted gas is assumed 
to be carried away with protostellar jets. The latter are 
triggered only after the formation of the central star. The 
fact that we use a sink cell means that our model cannot 
resolve the formation of binary (or multiple) stellar (or 
planetary) systems on spatial scales smaller than the size 
of the sink cell. 

2.1. Basic equations 

We make use of the thin-disk approximation to com- 
pute the gravitational collapse of rotating, gravitation- 
ally unstable cloud cores. This approximation is an ex- 
cellent means to calculate the evolution for many orbital 
periods and many model parameters and its justification 
is provided in Appendix JA) We note that in the thin- 
disk approximation all material from the envelope lands 
onto the outer disk regions. This is however a r eason- 
able assumption according to IVisser et al.l (J2009l ). who 
accurately calculated the gas trajectories in the infalling 
envelope and found that the bulk of the infalling material 
landed onto the disk's outer edge. 

The basic equations of mass, momentum, and energy 
transport in the thin-disk approximation are 



dt 



= -Vp • i^Vp) , 



(1) 



-{^Vp)^[\/-{J:vp^Vp)l = -\/pV^j:gp^{\/-U)p, 

(2) 

de 

— + Vp • (evp) = -V{Vp-Vp)-A^r^ ( Vi;)^^, : Upp^ , 

(3) 
where subscripts p and p' refers to the planar components 
(r, (j)) in polar coordinates, S is the mass surface density. 
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e is the internal energy per surface area, V = j_^ Pdz is 
the verticahy integrated form of the gas pressure P, Z is 
the radially and azimuthally varying vertical scale height 
determined in each computational cell using an assump- 
tion of local hydrostatic equilibrium (|Vorobvov & Basul 
l2Q09af ). Vp = VrV -\-V(f)(p is the velocity in the disk plane, 
9p = dr'f' + 9ci)^ is the gravitational acceleration in the 

disk plane, and Vp = rd/dr + <pr~^d/d(j) is the gradi- 
ent along the planar coordinates of the disk. The pla- 
nar components of the divergence of the stress tensor 
(V • n)p, symmetrized velocity gradient tensor (Vi;)^, 



IIpp/, and symmetric dyadic 



viscous heating (Vv) 
T^Vp (g) Vp are found according to the usual rules (see Ap- 
pendix C). 

The gravitational acceleration g^ includes the gravity 
of a central point-mass star (when formed), the gravity 
of an inner disk (r < rgc), and the self-gravity of a cir- 
cumstellar disk and envelope. The latter component is 
found by solving for the Poisson integral 



^(r,( 



-G 



prout 



n2TT 

Jo 



T.{r',(t)')d(t)' 



2rr' cosf 



where Tout is the radial position of the computational 
outer boundary, or, equivalent ly, is the initial radius of a 
cloud core. This integral is calculated using a FFT tech- 
nique which applies the 2D Fourier convolution the orem 
for polar coordinates (see lBinnev fc Tremainelfl987l Sect. 
2.8). 

2.2. Viscosity 

Viscosity in circumstellar disks may be an important 
source of mass and angular momentum transport and 
heat production. The best candidate to date is turbu- 
lent viscosity induced by the magneto-rotational instabil- 
ity (|Balbus fc Hawlevlll99lf ). though other mechanisms 
such as nonlinear hydrodynamic turbulence cannot be 
completely eliminated due to the large Reynolds num- 
bers involved. We make no specific assumptions about 
the source of turbulence and parameterize the magni- 
tude of kinematic viscosity using a modified form of the 
a-prescription 

V = ac^Z Tcx{r), (5) 

where Cg = 7P/S is the square of effective sound speed 
calculated at each time step from the model's known V 
and S. The function Toc{r) = 27r-^ td^n-^ [{r^/ry^] 
is a modification to t he usual a-prescription of 
IShakura fc SunvaevI (| 19731 ) that guarantees that the tur- 
bulent viscosity operates only in the disk and quickly 
reduces to zero beyond the disk radius r^- The latter is 
determined using a typical density for the disk to enve- 
lope transition, Sd2e = 0.1 g; cm~^, and the radial gas 
velocity (see .Vorobvovl 120101 for details). 

In this paper, we use a spatially and temporally uni- 
form a, with its value set to 0.005 in most models. This 
choice is based on our recent work (|Vorobvov fc Basul 
l2009a), wherein we have studied numerically the secular 
evolution of viscous and self-gravitating disks. We found 
that if circumstellar disks around solar-mass protostars 



could generate and sustain turbulence, then the tempo- 
rally and spatially averaged a should lie in the range 
10-^ - 10-^ Smaller values of a (< 10""^) have little 
effect on the resultant disk structure and mass accretion 
history, which, in this case, is totally controlled by disk 
gravity. Larger values {a > 10~^) destroy circumstel- 
lar disks during less than 1.0 Myr of evolution and are 
thus inconsistent with mean disk lifetimes of the order 
of 2-3 Myr. Nevertheless, a may vary in time and have 
greater values in the EPSF (the duration of which is usu- 
ally much shorter than 1 Myr). The effect of varying a 
is briefiy discussed in Section 15.31 

Viscosity enters the basic equations via the viscous 
stress tensor 11 expressed as 

U = 2^iy (\/v - i(V • v)e] , (6) 



where e is the unit tensor. We note that we take no 
simplifying assumptions about the form of 11 apart from 
those imposed by the adopted thin-disk approximation. 

2.3. Energy balance 

Equation (J3]) for the internal energy balance includes 
the usual compressional term V (Vp ■ Vp)^ radiative cool- 
ing A, heating due to stellar/background irradiation F, 
\^) and viscous heating {\/v)pp' : Upp' . We assume that 

the heat generated in the disk interior due to viscos- 
ity and shocks is transported to the disk surface by ra- 
diation, which escapes from the disk surface at a rate 
per unit area 2cfT^^. This means that we neglect other 
possible sources of heat transport such as convection. 
We then make use of the diffusion approximation and 
link the effective surface temperature Teff with the mid- 
plane temperature of gas T^p via the following relation 
TJj = 8Tj^p/(3r), where r is the optical depth (jHubenvl 
I1990I : iJohnson fc Gamnii3l2003[ ). Finally, we substitute 
T~^ with r/(l -h T^) to allow for a smooth transition 
between the optically thic k and optically thin regimes 
([Johnson fc Gammiell2003f ). The resulting cooling func- 
tion is described as 

A = J-c<tT^p-^, (7) 



where <j is the Stefan-Boltzmann constant and J^^ = 
2 + 20tan~-^(T)/(37r) is a function that secures a correct 
transition between the cooling function in the optically 
thick regime Athick = 16crT^p/(3r) and the optically 
thin one Athin = ^crTj^pT. We use frequency-integrated 
opacities of iBell fc Lid (|l994f ). which are smoothed at 
the principal opacity transitions to allow for iterative so- 
lution methods to converge quickly. 

Heating due to stellar and background irradiation is 
treated assuming that this process operates in the oppo- 
site direction to that of radiative cooling, i.e., radiation 
from the central star and natal molecular cloud hits the 
surface and diffuses down to the midplane where it trans- 
forms into heat. This allows us to express the heating 
function as 



F = JvCrT:^^ 



(8) 



•1 + t2' 

where Tjrr is the irradiation temperature at the disk sur- 
face determined by the stellar and background black- 
body irradiation as 

firr(r) 



-^ irr ^ b, 



(9) 
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Fig. 1. — Schematic representation of the numerical model. See 
the text for a detailed explanation. 

where Tbg is the uniform background temperature (in our 
model set to the initial temperature of the natal cloud 
core) and Firr(r) is the radiation flux (energy per unit 
time per unit surface area) absorbed by the disk surface 
at radial distance r from the central star. The latter 
quantity is calculated as 



i^irr(r) 



A ^* 
^irrT 2 ^^^^" 



(10) 



where L* is the stellar luminosity, 7irr is the incidence an- 
gle of radiation arriving at the disk surface at radial dis- 
tance r, and Ai^r is a time-dependent factor that accounts 
for the attenuation of stellar radiation in the EPSF (see 
Appendix [Bl for more details). 

The stellar luminosity L* is the sum of the accretion 
luminosity Laccr = GM*M/(2r*) arising from the grav- 
itational energy of accreted gas and the photospheric 
luminosity Lph due to gravitational compression and 
deuterium burning in the star interior. The stellar 
mass M* and accretion rate onto the star M are de- 
termined self-consistently during numerical simulations 
via the amount of gas passing through the sink cell. 
The stellar rad ius r^ is calculat ed usi n^ an approxima- 
tion formula of iPalla & Stahlerl (|1991[ ), modified to take 
i nto account the formation of the first molecular core 
(|Masunaga fc Inutsukall2QQQf ). More specifically, we as- 
sume that during 2 x 10^ yr after the formation of the 
central protostar, the stellar radius is r* = 5 AU. Then, 
the second atomic core forms and the stellar radius is 
determined as 

f 2.5 Ro M* < O.4M0, 

r* = <^ 2.5 + 4.2(M* - 0.4) R© OAMq < M* < LOM©, 
[ 5.0 Ro M* > 1.OM0. 

(11) 
Transition between these two modes is smoothed over a 

period of 0.5 x 10^ yr. 

The photospheric luminosity Lph is taken from the pre- 
main sequence track s for the low-mass star s and brown 
dwarfs calculated bv lD'Antona fc Mazitellil (|1997D . Un- 
fortunately, the stellar age in these tracks is difficult to 
relate with the actual physical evolution time in numer- 
ical simulations of gr avitational collaps e. A common 
practice starting from IM vers et al.l (|2005[ ) is to add ^offset 



to the times of the pre-main sequence tracks to account 
for the delay between the onset of cloud core collapse 
and the zero-time of these tracks. Indeed, after collapse 
begins, the forming star must wait for some time before 
the luminosity due to contraction and deuterium burn- 
ing (as described by D'Antona fc Mazitelli (1997)) will 
begin. The exact value of toffset is however uncertain be- 
cause it would certainly depend on the initial conditions 
in a cloud core such as the gas temperature, density en- 
hancement, strength of magnetic fields, etc. 

Fortunately, we accurately follow the pre-stellar col- 
lapse phase and can actually determine the time tfc that 
it takes for a cloud core to reach an optically thick den- 
sity of order 10^^ cm~^ in its interior and start forming 
the first (molecular) hydrostatic core. We then assume 
that the zero-time of D'Antona fc Mazitelli's tracks corre- 
sponds to the onset of the formation of the second atomic 
core tsc , which follows the f ormation of the first core afte r 
approximately 2 x 10^ yr (|Masunaga fc Inutsukal l2Q00f ) . 
With these assumptions in mind, the photospheric lumi- 
nosity L^^pYi is set to zero for t < ts.c. = tj.c. + 2 x 10^ yr 
and then is calculated according to lD'Antona fc Mazitellil 
(■1997. ) with the zero-time of their tracks corresponding 
to ts.c. in our numerical simulations. We note that the 
D'Antona fc Mazitelli's tracks do not cover the very 
early phases of stellar evolution. Therefore, we have 
used a power-law expression to extrapolate to times ear- 
lier than those included in the pre-main sequence tracks, 
Lph = ^ph,o(V^o)^5 where to is the earliest time in the 
tracks and Lph, is the pre-main sequence luminosity at 
this time. 

Viscous heating operates in the disk interior and is 
calculated using the standard expression {\/v)pp' : Upp' 
(see Appendix C). We note that we use the most general 
expression for viscous heating and take none of the pop- 
ular simplifying assumptions (such as disk axisymmetry) 
apart from those imposed by the thin-disk approxima- 
tion. Heating due to shock waves is taken into account 
via compressional heating V (Vp • Vp) and artificial vis- 
cosity. The latter is implemented i n the code u si ng th e 
standard prescription of Richtmeve r fc MortonI (|1957[ ). 
The gas pressure V and internal energy per surface area 
e are related via the ideal gas law 7^ = (7 — 1) e, with 
the ratio of specific heats 7 = 7/5. A more detailed 
appro ach would be to implement a variable 7 as in e.g. 
Forga net al.l (J2009l ). However, a rigorous realization of 
this mechanism requires calculating the excitation lev- 
els of main atomic and molecular coolants and is beyond 
the limits of the current paper. We explored the effect of 
varying 7 in our previous paper in th e context of poly- 
tropic disks (|Vorobvov fc Basul l2006f ) and showed that 
the 7 = 5/3 case was usually characterized by disks less 
prone to fragmentation. 

In Figure [1] we summarize our model by drawing a 
schematic picture of the main model ingredients in the 
EPSF. A central star is surrounded by a disk which ac- 
cretes matter from a collapsing natal cloud core. The 
infalling material lands onto the disk outer edge and is 
transported toward to the inner disk boundary by a com- 
bined action of gravitational and viscous torques. Mass 
accretion onto the star, along with stellar compression 
and deuterium burning, give rise to stellar irradiation, 
part of which is absorbed by the flaring disk surface and 
is transformed into heat in the disk interior. Another 
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TABLE 1 

Model parameters 



Model 


^ 


^0 


TQ 


Mci 


TbK 


a 


reference 


1.3 X 10-2 


2.7 


1540 


0.70 


10 


5 X 10-3 


lower- /3 


2.8 X 10-3 


1.2 


1640 


0.73 


10 


5 X 10-3 


Mci = 0.16 


1.3 X 10-2 


12 


340 


0.16 


10 


5 X 10-3 


Mci = 0.23 


1.3 X 10-2 


8 


514 


0.23 


10 


5 X 10-3 


Mci = 0.92 


1.3 X 10-2 


2 


2060 


0.92 


10 


5 X 10-3 


Tbg = 10 


1.3 X 10-2 


1.8 


2780 


1.2 


10 


5 X 10-3 


Tbg = 20 


1.3 X 10-2 


4.8 


1200 


1.1 


20 


5 X 10-3 


Tbg = 30 


1.2 X 10-2 


8.2 


860 


1.15 


30 


5 X 10-3 


a = 


1.3 X 10-2 


2.7 


1540 


0.70 


10 





a = 0.05 


1.3 X 10-2 


2.7 


1540 


0.70 


10 


5 X 10-2 



Note. — All masses are in Mq^ temperatures in Kelvin, dis- 
tances in AU, and angular velocities in km s— -"^ pc--*^. 

source of external heating is the background irradiation 
from the natal molecular cloud. The heat generated in 
the disk interior by viscosity and shocks is transported 
to the disk surface by radiation. The latter escapes from 
the disk surface giving rise to the only global cooling 
mechanism in our model. 

2.4. Solution procedure 

Equations ([I])-® are solved in polar coordinates (r, (j)) 
on a numerical grid with 512 x 512 grid zones. The radial 
points are logarithmically spaced. The innermost grid 
point is located at the position of the sink cell rgc = 
6 AU, and the size of the first adjacent cell varies in 
the 0.07-0.1 AU range depending on the cloud core size. 
This corresponds to a radial resolution Ar=l.l-1.6 AU 
at 100 AU. The outer boundary is reflecting. 

We use the method of finite differences with a time- 
explicit solut ion procedure similar i n methodology to the 
ZEUS code (IStone fc NormanlfT99 2). The advection is 
treated using the van Leer interpolation scheme. It is 
well known that cooling and heating time scales may be- 
come much shorter than the dynamical time scale, which 
would result in prohibitively small time steps. Therefore, 
the update of the internal energy per surface area e due 
to cooling A and heating T is done implicitly using the 
Newton- Raphson method of root finding, complemented 
by the bisection method where the Newton-Raphson it- 
erations fail to converge. The accuracy is guaranteed by 
not allowing e to change more than 30% over one time 
step. If this condition is violated in a particular cell, we 
employ subcycling for this cell, i.e., the solution is sought 
with a local time step that is smaller than the global time 
step by a factor of 2. The local time step may be further 
decreased until the desired accuracy is reached. 

The viscous force and heating terms in Equations (J2j) 
and (|3]) are implemented in the code using an explicit 
finite- difference scheme. This is found to be adequate for 
a < 0.01 because other terms (usually, the azimuthal ad- 
vection) dominate in the Courant condition that controls 
the time step. However, for higher values of a we find 
that the viscous terms start to impose strict time step 
limitations and an implicit scheme is desirable in order 
to extend numerical simulations to the Class II phase of 
stellar evolution. A small amount of artificial viscosity is 
added to the code to smooth out shocks. The associated 
artificial viscosity torques integrated over the disk area 
are negligible in comparison with gravitational torques. 
Occasionally, however, the shocks may become strong 



enough to impose strict limitations on the Courant con- 
dition, which results in a considerable decrease in the 
time step of integration. In this case, we use subcycling 
in the same manner as we do for the internal energy up- 
date due to cooling/heating. 

2.5. Initial conditions 

Initially, cloud cores have surface densities E and an- 
gular velocities Vt typical for a c ollapsing, a xisymmetric, 
magnetically supercritical core (iBasul[l997[ ): 



ro^o 



\/r'^ + rl 



n = 2n 



»( 



ro\^ 



-) 




(12) 



(13) 



where Qq is the central angular velocity and ro is the 
radius of central near-constant-density plateau defined as 
ro = vAc'^/{7rGTio)' We note that the above form of the 
column density is very similar to the integrated column 
density of a Bonnor-Ebert sphere (Dapp & Basu 2009|)- 
Furthermore, equation ([12]) at large radii r ^ ro leads to 
the gas volume density distribution p = Ac'^/{27rGr'^)^ if 
it is integrated in the vertical direction assuming a local 
vertical hydrostatic equilibrium, i.e., p = S/(2Z) and 
Z = c^/{7iGTi). This means that our initial gas surface 
density configuration can be considered to have a factor 
of A positive density enhancement compared to th at of 
the singular isothermal sphere psis = Cs/(27rGr^) (|Shul 
Il977j ). Throughout the paper, we use A = 1.2. 

Cloud cores are also characterized by the ratio of ro- 
tational to gravitational energy j3 = £^rot/|^grav|, where 
the rotational and gravitational energies are calculated 



as 
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(14) 
Here, a^ = ^ r is the centrifugal acceleration, and rout 
is the outer cloud core radius. The adopted values of /3 
lie within the limits inferred by ICaselli et al.l (|2002f ) for 
dense molecular cloud cores, /3 = (10~^ — 0.07). Cloud 
cores are initially isothermal, with the uniform gas tem- 
perature taking values between Tinit = 10 K and 30 K, 
depending on the model. In addition, every model core 
is characterized by a distinct ratio rout/^o = 6 in order 
to generate gravitationally unstable truncated cores of 
similar form. 

For the in-depth analysis, we consider a model with 
Mci = 0.7 Mo, /3 = 1.3 X 10-^ A = 1.2, a = 5 x 10"^ 
(the viscous a-parameter), and Tinit = 10 K. These and 
other model parameters are summarized in Table [H This 
model (hereafter, the reference model) is chosen solely 
because it best represents the main characteristics of disk 
fragmentation in the embedded phase of star formation. 
Other models will be introduced as the need arises. 

3. GRAVITATIONAL INSTABILITY AND DISK 
FRAGMENTATION 

Theoretical and numerical studies of the evolution of 
protostellar disks indicate that disk fragmentation is a 
complicated phenomenon, which can be influenced by 
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Fig. 2. — Gas surface density distribution (g cm ^, log units) in the reference model at several time instances after the formation of 
the central star (located in the coordinate center). Two bottom rows also shows a gas velocity field superimposed onto the surface density 
distribution. The vertical arrow in the bottom-right panel has a dimension of 5 km s"-*^. 



both the internal disk physics and external environ- 
ment. The latter may influence the disk susceptibil- 
ity to fragmentation directly (through, e.g., disk irra- 
diation) or indirectly by setting the initial conditions in 
cloud cores that favor or disfavor fragmentation in subse- 
quently formed disks. The following four criteria for disk 
fragmentation are best studied and their significance is 
well established. 

1. The ratio of the local cooling time tc = e/A to 
the local dynamical time Q~^ is smaller than a 



few, 



i.e., tc.n < C 



iMejia et al.l l2QQ5f ) . The actual value of C may vary 
depending on the physical conditions in the disk, 
e.g., C may depend on the disk thickness, chemi- 
cal composition, dust content, etc. In the following 
text, we will refer to the dimensionless quantity tc^ 
as the ^-parameter and adopt ^ = 1 as a fiducial 
critical value. 

2. The Toomre criterion Q = cS^/ijiGT^) for a Kep- 
lerian disk is smaller than some critical value Qcr, 
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usually taken to be unity ( jToomrd 1 1964 1). Here 
again, Qcr may depend on the physical conditions 
and may vary by a factor of unity. The Toomre 
criterion implies that the gas surface density S 
should be sufficiently high for a disk to fragment. 
This criterion, along with Q < 1-3, is often in- 
voked when analy zing the disk susceptibility to 
fragmentation (e.g. lRafikovll2QQ5T ). There is, how- 
ever, a catch — too high S may prevent fragmen- 
tation due to in creased opacity and cooling time 
(|Nero fc Bjorkm an 2009). In other words, there 
exists minimum and maximum values of S between 
which the instability and fragmentation are ex- 
pected to occur. This means that any numerical 
simulation that starts from a pre-defined star /disk 
system with some disk-to-star mass ratio may run 
the risk of not revealing disk fragmentation if the 
initial S is too high. This is generally not a prob- 
lem in numerical simulations that form disks self- 
consistently (such as our own), because during the 
disk formation phase E naturally increases from 
low toward higher values and the disk may pass 
through the unstable phase. 

3. The amount of rotation in the natal cloud 
core should be sufficiently large in order to 
form extended an d mas si ve protostella r disks 
(IVorob vov fc Basul l2006l: jEZatter et al.l 120081: 
VorobvQ^ l2009bl : iRice at al.l l2010l : iMachida et all 
201Q|). 
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4. The time-averaged rate of mass accretion onto the 
disk (Md) is greater than the time- averaged mass 
accretion onto the star (M*) so that S quickly 
increase s with time and may reach the unstable 
regime (lA^robvov fc Bas^ I2006L [2QQ3 ; IVorobvovl 
[2009bl : iKratter et al.l l2 1 Ol : iBd^ 120091 ). 

We analyze the significance of these four criteria for 
disk fragmentation using our reference model. Figure [2] 
shows a series of images of the gas surface density (in 
g cm~^, log units) in the inner 1000 AU at different times 
since the formation of the central star. The disk begins 
to form at t ^ 0.08 Myr and by t = 0.13 Myr a well- 
developed spiral pattern and several dense clumps are 
clearly visible. The clumps are almost always located 
in the spiral arms, suggesting that they form via frag- 
mentation of the densest and coldest arms. Most frag- 
ments, however, do not live long. They are driven into 
the disk inner regions and through the sink cell (and 
probably onto the star) but other fragments take their 
place. Some of them are massive enough to host mini- 
disks of their own. Typical fragment masses lie in a wide 
range from several Jovian masses to low- and interme- 
diate mass brown dwarfs. The mass spectrum of the 
fragments depends on the disk and cloud core properties 
and may vary from model to model. 

The disk slowly grows in mass and size due to mass 
loading from the envelope (most of which is off the spatial 
scale in Figure [2j). The disk structure is rather irregular, 
particularly in the early evolution. The gas velocity field 
in the bottom rows of Figure [2] reveals large non-circular 
motions, contractions, and expansions caused by ongoing 
angular momentum redistribution between the fragments 





y^^ 



J 



a 





(7 



i?5>^:: 



r 




-500 



500 -500 

Radial distance (AU) 



500 



Fig. 3. — Gas surface density distribution (left column, g cm~^, 
log units) and the spatial distribution of the ^-parameter (right 
column, log units) in the reference model at three typical times 
after the formation of the central star. Red contour lines comprise 
gravitationally unstable regions according to the Toomre criterion, 
Q < \. Fragmentation is supposed to occur in the regions where 
Q < 1 and Q < 1 simultaneously. 



and the rest of the disk (in particular, between the frag- 
ments and spiral arms). The disk in this early phase of 
evolution is most certainly not in a steady state and ap- 
proximating the early disk evolution using a steady-state 
concept may be misleading. 

Figure [2] reveals that the disk in the reference model is 
readily susceptible to fragmentation in the early evolu- 
tion. How does the model comply with the four fragmen- 
tation criteria outlined above? Are all four conditions 
satisfied? We start with examining the importance of 
criteria 1 and 2 and search for any disk regions that are 
simultaneously characterized by both Q < 1 and Q < 1. 
Figure [3] presents several typical gas surface density dis- 
tributions (left column, g cm~^, log units) and the spatial 
distribution of the corresponding ^-parameter (right col- 
umn, log units). In the latter case, those regions that cool 
sufficiently fast for fragmentation to take place (^ < 1) 
are plotted with blue, while slowly cooling regions with 
Q > 1 are plotted with red. Disk regions shown with 
white are near the border of stability, Q = 1. The black 
contour lines dehneate the regions of the disk that are 
prone to fragmentation according to the Toomre crite- 
rion, Q < 1. It is clearly seen that there are regions in 
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Fig. 4. — Same as in Figure [3] but for the lower-^ model, /5 = 
2.8 X 10""^. Note the lack of fragmentation. 

the disk where the first two criteria for fragmentation, 
G = tc^ < 1 and Q < 1, are fulfilled simultaneously. 
These are the fragments, especially those located in the 
outer disk regions, and certain parts of the spiral arms. 
It is seen that favorable sites for fragmentation lie prefer- 
ably at large radii, implying that many fragments form at 
r > 100 AU from the star but are driven later in the in- 
ner regions via exchange of angular momentum with the 
disk and, especially, with the spiral arms. This migration 
phenomenon"^ wa s demonstrated by us in the context of 
barotropic disks (|Vorobvov fc Basull2006f ). 

Figure [3] demonstrates that criteria 1 and 2 for disk 
fragmentation are fulfilled in the reference model. What 
about the other two criteria? Criterion 3 is essentially an 
initial condition imposed on the cloud core which states 
that the rate of cloud core rotation should be sufficiently 
high for disk fragmentation to take place. To investi- 
gate the importance of this condition, we consider an- 
other model that is similar to the reference model but 
has a smaller initial rotation rate (hereafter, lower- /3 
model). In particular, we set the ratio of rotational to 
gravitational energy to /3 = 2.8 x 10~^ (in contrast to 
P = 1.3 X 10~^ in the reference model) by decreasing 
the value of Qq in Equation ([13]). The resulting distri- 
butions of the gas surface density S, ^-parameter, and 

^ The animation of this migration process can be downloaded at 
www.ap.smu.ca/~vorobyov/ 



Toomre parameter Q are shown in Figure ID The layout 
of the figure is the same as that of Figure [3]but the spatial 
scale is different. It is evident that the lower- /3 model has 
no well-defined fragments, though some transient density 
enhancements within the spiral arms are visible. The 
lack of disk fragmentation is not surprising — there are 
hardly any regions in the disk where criteria 1 and 2 are 
satisfied simultaneously. In fact, by t = 0.3 Myr, the 
disk lacks regions with Q < 1 and regions with Q < 1 
are mostly located near the disk outer edge where in- 
tense cooling of the shocked gas (due to accretion from 
the envelope) takes place. 

There are two major factors that work against disk 
fragmentation in the lower- /3 model. First, the disk size 
is considerably smaller than that of the reference model 
due to a smaller centrifugal radius r^ = Vfir^ /GM{r). 
Smaller disks are subject to a stronger stabilizing influ- 
ence of stellar irradiation. Second, the disk mass in the 
lower- /3 model is on average 20% smaller than that of the 
reference model, which also increases the disk stability 
against fragmentation in the lower-/3 model by raising 
the value of Q. In addition, small disks may be opti- 
cally thick and thus cooling too slow to fragment (e.g. 
iRice at al.lf2009l : IClarkell2009f ). 

The above analysis indicates that the initial conditions 
in a natal cloud core (in particular, the amount of rota- 
tion), are of considerable importance for the future disk 
evolution. In models with low /3, the resulting disks are 
unlikely to fragment due to small disk sizes and masses. 
In this sense, criterion 3 is a necessary condition for disk 
fragmentation but not a sufficient one. As will be demon- 
strated later, disk propensity to fragment also depends 
on other factors such as magnetic fields, initial cloud core 
temperature and mass, etc. In this context, it is difficult 
to provide reliable estimates as to the exact amount of 
rotational energy (as specified, for example, by the ra- 
tio P of the rotational to gravitational energy) that a 
cloud core needs in order to produce disks capable for 
fragmentation. Therefore, we believe that providing any 
critical values of /3 for disk fragmentation may be mis- 
leading unless exact initial conditions in cloud cores are 
specified. 

In the following section, we will consider mass accre- 
tion rates onto the disk and the star and discuss the 
significance of criterion 4 for disk fragmentation. 

4. THE BURST MODE OF ACCRETION 

Fragments that form in the disk pass through the sink 
cell as they migrate into the inner disk via exchange of 
angular momentum with the spiral arms. The ultimate 
fate of these fragments is uncertain and is largely de- 
pendent on how quickly they can contract from their 
initial size of several AU to a planetary size to avoid 
tidal destruction. The contraction time for a Jupiter- 
mass clump to reach a central temperature of 2000 K, 
i.e., the temperature required to dissociate H2 to trig- 
ger rapid collapse , may be as long as a few x 10^ yr 
(|Helled et all 12006V Considering a fast timescale of in- 
ward radial migration in the embedded phase — a few 
tens of orbital periods — we believe that most of these 
fragments^ will be tidally destroyed when approaching 

^ The most massive fragments may survive and form giant plan- 
ets or brown dwarfs on close orbits. 
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Fig. 5. — Mass accretion rate onto the star (top) and stellar 
luminosity (bottom) as a function of time elapsed since the begin- 
ning of collapse in the reference model. In particular, the solid and 
dashed lines in the bottom panel show the accretion and photo- 
spheric luminosities, respectively. 



the central star, thus converting its gravitational energy 
to the accretion luminosity and producing an FU-Ori-like 
luminosity burst. This phenomenon is called the burst 
mode of accretion and it has been extensively studied 
by us for the case of barotropic disks (Vorobyov & Basu 
l2QQ5Ll2QQ6l ). Here, we confirm that a more accurate treat- 
ment of disk thermodynamics does not qualitatively af- 
fect our earlier conclusions. However, as our recent simu- 
lations of barotropic disks have shown, some of the frag- 
ments that form in the late embedded phase may survive 
and evolve eventually into giant planets on wide orbits 
(|Vorobvov fc Basul HoiQ). 

The instantaneous mass accretion rate from the disk 
onto the star M* is found in our model as the mass pass- 
ing through the sink cell per one time step of integra- 
tion (which in physical units is usually equal to 10-20 
days). We also calculate the instantaneous mass accre- 
tion rate onto the disk from the infalling envelope M^ 
as the mass passing (per one time step of integration) 
through a radial annulus located just outside the disk 
outer edge. Figure [5] presents the time evolution of the 
mass accretion rates and luminosities in the reference 
model. In particular, the top panel shows M*, while the 
bottom panel — accretion luminosity Laccr (solid line) and 
photospheric luminosity Lph (dashed line). 

In the pre-stellar phase, M* is negligible but quickly 
rises to ~ 10 ~^ Mq yr~^ when the gas volume density 
in the sink cell exceeds 10^^ cm~^ and a central stellar 
core begins to form at t ~ 0.08 Myr after the onset of 
collapse. The subsequent short period of near-constant 
accretion corresponds to the phase when the infalling en- 



velope lands directly onto the forming star^. A sharp 
drop in M* follows shortly, which manifests the beginning 
of the disk formation phase. Subsequently, the infalling 
envelope accretes onto the forming disk rather than di- 
rectly onto the star. This transient drop in M occurs 
due to the fact that the disk mass is initially too small 
to drive a substantial accretion rate onto the star either 
due to viscous or gravitational torques. As the evolu- 
tion proceeds, the disk accretes mass from the infalling 
envelope and a qualitatively new phase of mass accre- 
tion ensues, in which M* shows variability by several 
orders of magnitude. Short episodes of high-rate accre- 
tion (caused by the passage of disk fragments through 
the sink cell) are followed by longer periods of low-rate 
accretion (caused by a temporary disk expansion and sta- 
bilization). This highly variable accretion makes the star 
sporadically increase its total luminosity, as illustrated in 
the bottom panel of Figure [5l Several clear-cut luminos- 
ity outbursts with Laccr as high as 100 Lq and many 
more weaker bursts (solid line) are evident against the 
background of a near-constant photospheric luminosity 
with Lph ~ 1.0 Lq (dashed line). The stronger bursts 
may represent FU Orionis-like eruptions, typical for the 
early evolution of a protostar, while weaker ones may 
manifest EX Lupi-like eruptions (EXors), typical for the 
later evolution. We note that the exact time for the 
onset of the photospheric luminosity is rather uncertain 
and may shift to later times (see discussion in Section [6]), 
which would result in the early luminosity bursts being 
considerably stronger in amplitude. 

We can now verify if our reference model complies with 
criterion 4 for disk fragmentation outlined in the pre- 
vious section. This criterion requires that the rate of 
mass accretion onto the disk M^ be on average greater 
than that onto the star M*. Figure [6] presents the time- 
averaged mass accretion rates onto the star (M*) (solid 
line) and onto the disk (Md) (dashed line) as a function 
of time since the beginning of collapse. The averaging is 
done over a period of 15000 yr. In the early evolution 
{t < 0.2 Myr), (Md) is systematically greater than (M*) 
and this phase is characterized by the strongest burst ac- 
tivity. In the subsequent time period between 0.2 Myr 
and 0.3 Myr, both time- averaged accretion rates are of 
similar magnitude and the burst phenomenon persists, 
though with somewhat lesser frequency and amplitude. 
After t = 0.4 Myr, (Md) becomes systematically lower 
than (M*) and the burst activity in this late phase di- 
minishes. However, some small variations in M* persist 
even to later times. 

Let us define the end of the embedded phase and the 
onset of the Class II phase of star formation as the time 
when the envelope empties, and its mass Menv drops be- 
low 5-10% of the initial cloud core mass Md. The vertical 
lines in Figure [6] correspond to the evolution times when 
Menv/Mci = 0.1 (left) and Menv/Md = 0.05 (right). It is 
seen that (Md) > (M*) in the Class and I phases, while 
(Md) < (M*) in the Class II phase. Hence, disk fragmen- 
tation and the associated burst phenomenon are likely to 
take place in the embedded phase of star formation^ but 

^ In fact, this period is expected to be even shorter since Tsc ^ 
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Fig. 6. — Time- averaged accretion rates onto the star (solid 
line) and onto the disk (dashed line) versus time since the onset 
of collapse in the reference model. The vertical dotted lines mark 
the onset of the Class II phase as inferred from the ratio of the 
envelope mass to the initial cloud core mass, 0.1 and 0.05 for the 
left and right lines, respectively 

are unlikely later in the evolution simply because mass 
loading from the envelope diminishes in this phase. 

Figure [6] demonstrates that criterion 4 is fulfilled in 
the reference model. Is this criterion sufficient for disk 
fragmentation to take place? In Figure [7] we present the 
instantaneous accretion rates M* (top panel) and time- 
averaged accretion rates (bottom panel) in the lower-/3 
model introduced in the previous section. This model has 
a (roughly) five times smaller value of /3 = 2.8 x 10~^ 
as compared to that of the reference model and shows 
hardly any signs of disk fragmentation (see Figure H]). 
The lack of disk fragmentation manifests itself by a con- 
siderably weaker accretion variability than in the refer- 
ence model — there are only order-of-magnitude fiicker- 
ing in M* and one moderate accretion burst. However, 
when we turn to the time-averaged accretion rates (bot- 
tom panel), we see that (Menv) (dashed line) is actually 
greater than (M*) (solid line) in the EPSF, indicating 
that criterion 4 for disk fragmentation is fulfilled in the 
lower- /3 model. This example convincingly demonstrates 
that the fulfillment of criterion 4 is necessary but not 
sufficient for disk fragmentation to occur. The disk mass 
and radius in the lower-/3 model seem to be too small 
even in the case of a strong mass loading from the enve- 
lope. 



5. 



THE EFFECT OF INITIAL CONDITIONS ON THE 
BURST MODE OF ACCRETION 



In Section [31 we have already demonstrated the impor- 
tance of rotation for the development of the burst mode 
of accretion in the early phases of stellar evolution. In 
this section, we study the effect that other initial condi- 
tions in collapsing cloud cores (such as cloud core mass 
and temperature, magnetic fields, etc.) may have on the 
strength and frequency of the bursts. 

5.1. Initial cloud core mass 

There is at least one good reason to believe that the 
initial mass of a cloud core should have a significant effect 
on the subsequent disk evolution — more massive cloud 
cores are expected to form more massive disks. This 
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Fig. 7. — Instantaneous mass accretion rate (top) and time- 
averaged mass accretion rates (bottom) in the lower-/3 model, 
which is characterized by a five times lower /5 = 2.8 x 10~^ than 
that of the reference model. In particular, the solid and dashed 
lines in the bottom panel show the time-averaged accretion rates 
onto the star and onto the disk, respectively. Vertical dotted lines 
mark the onset of the Class I (left) and Class II (right) phases. 

is simply because more massive cloud cores have larger 
sizes'" and, as a consequence, larger centrifugal radii rcf 
for any reasonable radial mass distribution. Hence, we 
can expect disks formed from more massive cloud cores to 
have a higher tendency for fragmentation and a stronger 
accretion variability. This effect h as been c onfirmed in 
the context of barotropic disks (Vo robvov fc B asu 2QQ6I : 
Vorobvov 2QQ91J). A si milar tendency was demonstrated 
by Kratter et al'1(|2QQ8l ). who showed that stars of greater 
mass tend to have disks that are more susceptible to 
fragmentation. 

Figure [8] presents the mass accretion rates onto the 
star (left column) and accretion and photospheric lu- 
minosities (right column) in three models with Mc\ = 
0.16 Mq (top row), Mci = 0.23 Mq (middle row), 
and Mci = 0.92 Mq (bottom row). In the following 
text, we refer to these models as the Mc\ = 0.16 Mq 
model, Mci = 0.23 Mq model, and M^i = 0.92 Mq 
model, respectively. Other parameters of these models 
are identical to the parameters of the reference model 
and are summarized in Table [H It is seen that mod- 
els with lower Md are characterized by a lower accre- 
tion variability, suggesting that the disk propensity to 
fragment declines with decreasing cloud core mass. The 
Mci = 0.16 Mq model exhibits hardly any (or very weak) 
accretion and luminosity bursts, with the photospheric 
luminosity dominating the total radiation fiux for most of 
the evolution. The burst mode becomes prominent in the 
Mci = 0.23 Mq model, which shows three well-defined lu- 

^ A cloud core may also increase its mass via density enhance- 
ment. 
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Fig. 8. — Mass accretion rate onto the star (left column) and 
stellar luminosity (right column) as a function of time since the 
beginning of collapse in the M^ = 0.16 Mq model (top), M^ = 
0.23 Mq model (middle), and M^ = 0.92 Mq model (bottom). In 
particular, the solid and dashed lines in the bottom panel show the 
accretion and photospheric luminosities, respectively. 

minosity outbursts. As the cloud core mass continues to 
increase, the burst frequency and intensity also increase 
and the Mc\ = 0.92 Mq model demonstrates multiple lu- 
minosity outbursts with Laccr ^ 10-100 Lq and several 
ones with Laccr > 100 L©, indicating the onset of vigor- 
ous gravitational instability and disk fragmentation. 

We point out that all three models have the same 
value of /3 = 1.3 X 10"^ yet the M^i = 0.16 Mq and 
Mci = 0.23 Mq models have a considerably weaker burst 
activity than the M^ = 0.92 Mq model. This exam- 
ple demonstrates the importance of the initial cloud core 
mass for the development of disk fragmentation and as- 
sociated burst mode of accretion. For disk fragmentation 
to take place, it is not sufficient for a cloud core to have a 
high initial rate of rotation — the initial cloud core mass 
should also be sufficiently high. We also note that as 
Mci increases in Figure [H the resulting total luminosity 
also increases but this does not suppress disk fragmen- 
tation. The growing disk mass outweighs the stabiliz- 
ing influence of stellar irradiation, at least for stars with 
M^ S, 1.0 M(7). Our conclusion is in line with that of 
iRice at ahl (|2010f ) who argue that the primary require- 
ment for disk fragmentation is large enough /3 to produce 
disks with radii large enough for fragmentation. Indeed, 
we may form disk of greater size not only by increasing 
/3 but also by taking a larger (and hence more massive) 
cloud core. 

5.2. Higher initial cloud core temperature 

In the reference model, we set the initial cloud core 
temperature to Tinit = 10 K. According to our model 
assumptions, this value is physically determined by the 
temperature of the background blackbody radiation Tbg, 
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Fig. 9. — Mass accretion rates (left column) and stellar lumi- 
nosity (right column) as a function of time since the beginning of 
collapse in the Tbg=10 K model (top), Tbg=20 K model (middle), 
and Tbg =30 K model (bottom). In particular, the black solid and 
red dashed lines in the left column present the instantaneous mass 
accretion rate onto the star and time-averaged mass accretion rate 
onto the disk, respectively, while these lines in the right column 
show the accretion and photospheric luminosities, respectively. A 
color version of this figure is available in the online journal. 



i.e., Tinit = ^bg- However, Tbg may be higher and this 
may influence the subsequent evolution of the cloud core 
in at least three ways. First, the rate of mass accretion 
onto the disk will be greater because Md is proportional 
to the cube of the sound speed. This effect will assist 
disk fragmentation. Second, the background radiation 
flux will grow and moderate the disk tendency to frag- 
ment by systeni atically increasing the disk temperature 
(jCai et al.ll2008f ). And lastly, an increased rate of mass 
accretion onto the disk may eventually lead to an in- 
creased rate of mass accretion onto the star, thus rais- 
ing the accretion luminosity and contributing to another 
factor against disk fragmentation. It is unclear a priori 
which of the three key factors would dominate the disk 
evolution. 

To study the effect of varying background temperature, 
we consider three models that have similar cloud core 
masses and rotation rates but different background tem- 
peratures: Tbg =10 K, Tbg=20 K, and Tbg=30 K. In the 
following text, we refer to these models as the Tbg = 10 K 
model, Tbg = 20 K model, and Tbg = 30 K model, re- 
spectively, and their parameters are listed in Table [H 
We specifically choose models with similar Mc\ and P in 
order to avoid any possible interference with the effects 
based on different cloud core masses and rotation rates 
considered in Section |4] and 15.11 respectively. Figure [9] 
presents mass accretion rates (left column) and luminosi- 
ties (right column) as a function of time since the onset of 
gravitational collapse in the Tbg=10 K model (top row), 
Tbg=20 K model (middle row) and Tbg=30 K model (bot- 
tom panel). In particular, the black solid and red dashed 
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Fig. 10. — Mass accretion rates onto the star (left column) and 
stellar luminosity (right column) as a function of time since the 
beginning of collapse in the a = model (top) and a = 0.05 
model (bottom). In particular, the solid and dashed lines in the 
right column present the accretion and photospheric luminosities, 
respectively. 

line in the left column are the instantaneous mass accre- 
tion rate onto the star M* and the time-averaged (over 
15000 yr) mass accretion rate onto the disk (Md), respec- 
tively. The black solid and red dashed lines in the right 
column are the accretion and photospheric luminosities, 
respectively. 

A comparison of the three models reveals that the 
Tbg=20 K model exhibits a vigorous burst activity 
comparable in strength and frequency to that of the 
Tbg=10 K model. However, the duration of the burst 
phase appears to be shorter in the higher- Tbg model. 
As we further increase the background temperature to 
Tbg =30 K, the burst activity decreases notably, yet there 
are two well-defined accretion and luminosity outbursts 
that reveal the disk is still prone to fragmentation. In 
fact, the magnitude of these bursts is much stronger than 
in the lower- Tbg models, indicating that a higher back- 
ground radiation favors the formation of more massive 
fragments (though in a much smaller quantity). As was 
expected from theoretical grounds, the photometric lu- 
minosity T*^ph is greater in models with higher Tbg, but 
so is the mass accretion rate onto the disk (Md) (at least 
in the early phase). It appears that an elevated mass 
accretion rate onto the disk outweighs the stabilizing in- 
fiuence of the background and stellar irradiation. We 
conclude that protostellar disks illuminated by the back- 
ground irradiation with temperatures of the order of 30 K 
(and probably higher) are still prone to fragmentation 
and development of the burst mode of accretion. 

5.3. The effect of viscosity 

The effect of disk viscosity on the burst mode of ac- 
cretion was_studi^d__b}^sin__t he context of barotropic 
disks (|Vorobyov fc Basul l2009a[). For the usual a- 
parameterization of IShakura fc SunvaevI (|l973f ) and tem- 
porally and spatially constant a, disks with lower values 
of a are disposed to stronger fragmentation and demon- 
strate a stronger burst mode of accretion. In addition, 
the accretion variability al so increases alo ng!; the line o f 
decreasing a (JVorobvov fc Basu.2009a, : ,Vorobvovll2009bl ) . 

In all models considered so far, we have adopted a = 



0.005. To see how different values of a could affect our 
conclusions, we run two models with a = and a = 0.05 
but other parameters identical to those of the reference 
model (see Table [T]). Figure [10] presents the mass ac- 
cretion rates onto the star (left column) and luminosi- 
ties (right column) in the a = model (top row) and 
a = 0.05 model (bottom row). As expected, the a = 
model demonstrates a vigorous burst activity, while the 
a = 0.05 model shows only one strong luminosity out- 
burst with Taccr in excess of 100 T©, with other outbursts 
characterized by Taccr ^ 20 Lq . We confirm that a factor 
of 10 increase in a does not suppress disk fragmentation 
completely. However, an additional strong source of mass 
transport via viscous torques reduces the disk mass and 
this acts to moderate the disk propensity to fragment. 

The problem with the a = 0.05 model is that 
it demonstrates hardly any accretion episodes with 
M* < 10"^ Mq yr"^ in the early 0.2 Myr of evo- 
lution. The lack of low-rate accretion in disks with 
g > 0.01 was also found in the contex t of barotropic disks 
(|Vorobvov fc Basull2009al : IVorQbvovll2009b) and th is con- 
fronts recent observations of lEnoch et aP (|2009f ). who 
find that a considerable fraction of Class I sources in 
young star-forming regions have inferred accretion rates 
below 10"^ Mq yr-^ We note that both the a = 
and a = 0.005 models show plenty of such low-accretion 
episodes. We therefore argue that disk viscosity in 
the embedded phase is unlikely to be characterized by 
a > 0.01. 

In a broader context of viscous {a < 0.01) versus non- 
viscous {a = 0) models, the former seem to yield accre- 
tion rates in the Class II (or T Tauri) phase that are a 
factor 2-3 greater than those of the non- viscous model 
(|Vorobvov fc Basul l2007l [2008f ). As a result, an addition 
of Qf-transport helps to bring Clas s II disk niasses in bet- 
ter agreement with observations (jVorobvovl l2009a[ ) . On 
the other hand, the early disk evolution (Class and I 
phases) is weakly affected by a-viscosity because mass 
and angular momentum transport in this st age is largely 
dominated by gravitational torques (V orobvov fc Basul 
2009a). This can also be seen from the comparison of 
Figure [5] with the top panels of Figure [TOl — there is little 
qualitative difference in the mass accretion history be- 
tween the a = and a = 0.005 models in the EPSF. 

6. MODEL LIMITATIONS 

In this section, we discuss several assumptions in our 
model that can potentially infiuence our results. 1) The 
onset of photospheric luminosity. As was discussed in 
Section 12. 3| the stellar age in D'Antona fc Mazitelli's 
(1997) pre-main sequence evolution tracks is difficult to 
relate to the physical evolution time in numerical sim- 
ulation of cloud core collapse. We have equated their 
zero-point time to the time when the second atomic core 
presumably starts to form in our numerical simulations. 
This may be a conservative assumption. In D'Antona 
fc Mazitelli models, the evolution generally begins from 
a central temperature of log(Tc) = 5.7, i.e., at a time 
instance just preceding deuterium burning, and it may 
take some time for the forming second core to ignite deu- 
terium burning in its interior. Hence, the photospheric 
luminosity may turn on somewhat later than assumed in 
our numerical simulations and this could actually act to 
increase the disk susceptibility to fragmentation. 
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Fig. 11. — Aspect ratio A of the disk vertical scale height to 
radius (Z/r) as a function of radius (r). The thick solid line shows 
A in the reference model at t = 0.19 Myr after the formation of the 
central star. The dashed line presents the aspect ratio derived from 
the following expression Z/r = 0.1(r/100 AU)^-^^, as sugge sted by 
detailed disk vertical structure modeling bv ID Alessio et all {1999). 

2) Accretion rate onto the star. In our models, the size 
of the sink ceh rsc=6 AU is larger than the stellar radius. 
The inner disk at r < 6 AU may add additional variabil- 
ity to the accre tion rates, in par ticular due to the ther- 
mal instab ility (iBell fc Lid Il994l| or magne to-rotational 
instability (fArmitage et alJl2QQll : IZhu et alJl2QQ9l ). These 
effects may somewhat alter the temporal behavior of the 
actual accretion rates onto the stellar surface and affect 
the accretion luminosity. However, our accretion rates 
are in good accord with those inferred for nearby star- 
forming regions and we believe that the actual accretion 
rates onto the stellar surface are not substantially differ- 
ent from those calculated in our modeling. 

3) Jet efficiency. Protostellar jets may evacuate a sub- 
stantial fraction of the accreting material, reducing the 
effective mass of the star as compared to the case with- 
out jets. In our modeling, we set the jet efficiency to 
10%, which means that the stellar mass is systematically 
lower by 10% than that of the non-jet case. This in fact 
promotes gravitational instability in the disk as the disk- 
to-star mass ratio is increased accordingly. However, the 
jet efficiency may be higher and amount to 30% and pos- 
sibly more (e.g., Shu et al. 1999). This would act to fur- 
ther destabilize the disk. 

4) Stellar wobbling. The position of the central star in 
our models is fixed in the coordinate center. However, 
the star may move in response to the non-axisymmetric 
gravitational field of the disk. Semi-analytic consider- 
ations suggest that this stellar wo bbling may a mplify 
gravitational instability in the disk (|Adams et al.l fl989). 
though the recent numerical hydrodynamics simula tions 
find this effect insignificant (Scott fc Durisenll2010f ). To 
implement such a mechanism in our models is however 
not easy due to the presence of singularity in the coor- 
dinate center on the polar grid. We plan to explore the 
effect of stellar wobbhng in a future study. 

5) Binary or multiple system formation. In addition to 
the clump formation that we see in the present models, 
we have also seen the formation of a binary companion 
(or multiple companions) in the outer disk in models with 
Mci > 1.7 Mq and /3 > 2.0 x 10"^ These models wih 



be presented in a future paper. We note that to fully 
capture binary formation in the outer regions with our 
logarithmic grid, we need even higher numerical resolu- 
tion than in the present study. 

5) Magnetic fields. Frozen-in magnetic fields moderate 
the burst a ctivity due to an effecti ve increase in the Q- 
parameter (jVorobvov fc Basull2006l ). For a spatially and 
temporally uniform mass-to-fiux ratio, the magnetic ten- 
sion acts as a simple dilution of gravity, thus effectively 
lowering the disk surface density, and the magnetic pres- 
sure is a multiple of the gas pressure, thus providing an 
effective increase to the disk sound speed. A more com- 
prehensive study of the effect of magnetic fields, includ- 
ing ambipolar diffusion and magnetic braking, is planned 
for a future paper. 

7. CONCLUSIONS 

We have revisit ed our original r esults on the burst 
mode of accretion ( Vorobvov fc Basu 2005, 20^|), paying 
special attention to the thermal processes in protostellar 
disks around low- mass protostars. Our new model takes 
into account radiative cooling from the disk surface, vis- 
cous and shock heating, and also stellar and background 
irradiation. Thanks to the use of the thin-disk approx- 
imation, we can run uninterrupted numerical hydrody- 
namics simulations from the prestellar phase to the early 
T Tauri phase, fully capturing the embedded phase of 
star formation (EPSF). We find the following. 

• The EPSF is likely the only episode of disk evolu- 
tion when disk fragmentation can take place. How- 
ever, disk susceptibility to fragmentation in this 
phase depends crucially on the initial conditions in 
a natal cloud core. 

• Higher initial core angular momentum and mass 
lead to the formation of more massive and ex- 
tended disks and, therefore, favor disk fragmenta- 
tion. On the other hand, a higher temperature of 
the background irradiation Tbg may moderate the 
disk propensity to fragment. In particular, higher 
Tbg appears to favor the formation of more massive 
fragments though in much fewer numbers. 

• A higher rate of mass infall onto the disk than that 
onto the star in the EPSF does not guarantee disk 
fragmentation if the disk is not sufficiently large 
and massive. 

• For disk fragmentation to occur, both the Toomre 
Q-parameter and ^-parameter (ratio of the local 
cooling time to the dynamical time) must be be- 
low some critical value, taken to be unity in this 
paper, confirming many previous studies on disk 
instability and fragmentation. 

• Most (but possibly not all) fragments that from in 
the EPSF are driven into the inner disk regions 
and probably onto the star, triggering mass ac- 
cretion and luminosity bursts similar in magnitude 
to those of the FU-Orionis-type and EX-Lupi-like 
stars. This burst mode of accretion is a robust 
phenomenon that is expected to exist in a variety 
of environments and for a variety of systems with 
different physical properties. The intensity of the 
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burst mode correlates with the disk propensity to 
fragment. 

• Fragmenting disks drive highly variable accretion 
rates onto the star ranging from 10~^ Mq yr~^ to 
10~^ Mq yr~^. Protostellar disks that are gravi- 
tationally unstable but stable to fragmentation are 
characterized by a considerably weaker accretion 
variability with only an order of magnitude flicker- 



The intensity of the burst mode of accretion is sen- 
sitive to the amount of a- viscosity present in proto- 
stellar disks and appears to subside with increasing 
a. The lack of strong variability in disks with a spa- 
tially and temporally uniform a > 0.05 contradicts 



observations (e.g. lEnoch et al.l l2009f ) and renders 
such disk not viable. 
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APPENDIX 
A. THE THIN-DISK APPROXIMATION 

The thin-disk approximation is well justified as long as the aspect ratio A = Z/r of the disk vertical scale height Z 
to radius r does not considerably exceed 0.1. In a Keplerian disk, Z = Cs/^ and noticing that the angular velocity is 
Q = {GM^/r^y^'^ and the sound speed is Cg < QcrT^GTi/ft^ the aspect ratio can be expressed as 



A< 



■M^{r) 



CM, 



(Al) 



where M(^{r) = j T^{r^(j))r dr d(j) is the disk mass contained within radius r, M, is the mass of the central star, Qcr 
is the critical Toomre parameter, and C is a constant, the actual value of which depends on the gas surface density 
distribution S in the disk. For a disk of constant surface density, C is equal unity. However, circumstellar disks are 
characterized by surface density profiles declining with radius. For the scaling E ex r~^-^ typical for our disks, C = 4. 
Adopting further Qcr = 2 and M^{r)/M, = 0.5, which are typical tapper limits in our numerical simulations, we obtain 
A < 0.25. This analysis demonstrates that the thin-disk approximation is certainly valid in the inner regions where 
Md(r)/M* is small, but may become only marginally valid at large r where Md(r)/M* would approach its maximum 
value. 

The azimuthally- averaged radial distribution of the aspect ratio A = Z/r in the reference model at t = 0.19 Myr 
after the formation of the central star in shown by the solid line Figure [Til The vertical scale heig!:ht Z is calculate d 
assuming a local vertical hydrostatic equilibrium in the disk using the method described in lVorobvov fc Ba"sul (|2009aj ). 
Figure [TT] reinforces our analytical estimates and demonstrates the thin-disk approximation is certainly obeyed in the 
disk. Our disks rarely exceed 1000 AU in radius and the corresponding aspect ratio is kept in the 0.1-0.4 limits. Only 
at radial distances well in excess of 1000 AU may the thin-disk approximation be violated. 

B. RADIATION FLUX FROM THE CENTRAL STAR 

In order to calculate the radiation fiux from the central star Fi^r at a given radial distance r using Eq. (|9]), one needs 
to know the incidence angle of radiation arriving at the disk surface 7irr (i.e., the angle between the light rays and the 
perpendicular to the disk surface). For a fiaring disk, the cosine of 7irr can be expressed as 



cos 7ir 



COS airr COS /^irr (taU ai^r — taU /3irr) , 



(Bl) 



where cosofirr = dr/{dr'^ -\- dZ^^^^ •> cos/3irr = tI(t^ + Z^)-*^/^, tanairr = dZ/dr^ and tan/3irr = Z/r. In most cases, 
cosQfirr ~ 1 and cos/3irr ~ 1, siucc Z/v <C 1 (thin disk) and dZ/dr <C 1 (weak fiaring). Nevertheless, we use the 
complete expression for cos7irr. 

In the reality, the disk surface may not always be of the concave shape, so that tanairr — tan/3irr > and F\^^ > 0. 
If the shape of the disk surface becomes convex, i.e., tanairr — tan/3irr < 0, the irradiation fiux Firr becomes negative. 
Physically, this corresponds to the situation when part of the disk surface is shielded from the incoming radiation 
from the central star, for instance, by a local puffing of the disk. This may be potentially an important phenomenon. 
However, taking this effect into account self-consistently may require the use of full radiation transfer using ray tracing 
and is out of scope of the present paper. Therefore, to avoid this complicat ion, we make use of the detailed vertical 
structure models of irradiated accretion disks around T Tauri stars bv .D'Alessio et all (|l999f ). From their figure 1(b) 
(dashed curve) we have derived the following expression Z/r = 0.1(r/100 AU)^'^^ for the aspect ratio Z/r as a function 
of radial distance r, where 0.1 is the ratio Z/r at r = 100 AU and exponent 0.25 determines the degree of disk fiaring 
(for positive/negative exponents, the disk surface is concave/convex, respectively). We adopt this relation with a 
modification according to our model, i.e., we actually calculate the aspect ratio Z/r at 100 AU using the azimuthally 
averaged value of the vertical scale height Z. This would allow us to dynamically adjust the aspect ratio Z/r according 
to the actual disk thickness but keep the disk shape concave throughout the simulation. 
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Another effect that has to be taken into account is the attenuation of stehar irradiation by the infahing envelope 
in the embedded phase of star formation (EPSF). This is done by introducing a factor Ai^r in Eq. (|9]) calculated as 
Airr = Mci/(Menv + ^ci), whcrc Md is the cloud core mass (stays fixed) and Menv is the gradually decreasing envelope 
mass. In the early EPSF, Md ~ Menv and Ai^ ^ 0.5, while in the late EPSF, Menv -^ and Arr -^ 1. 

C. SUPPLEMENTARY MATHEMATICAL FORMULA 

For the convenience of the reader and for completeness, we provide the actual expressions for (V • n)^, (Vi;)pp, and 
[V • {TiVp Vp)] used in our paper. The components of V • 11 in polar coordinates (r, </>) are 



(v-n), 
(v-n). 



ld_ 
r dr 
d 



rUr 



1 d 



dr 



n. 



ld_ 

r d<p 



rllrd 



n. 



n. 



n 



r0 



(Cl) 
(C2) 



where we have neglected the contribution from off-diagonal components 11^;^ and Hcf^z • The components of the viscous 
stress tensor 11 in polar coordinates (r, (j)) can be found from Eq. ^ according to the usual rules. 

When calculating the symmetrized velocity gradient tensor Vv, only the following planar components are assumed 
to be non-zero: 



The symmetric dyadic T^Vp (g) Vp is a rank-two tensor expressed in polar coordinates (r, (j)) as 
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TiVrVr TiVrVcf) 

The planar components of [V • {T^Vp Vp)] can then be found using Eqs. (|C1|) and (|C2[) with 11 substituted by 



T.V 



p yy ^p- 

Finally, the viscous heating term (Vv) , : 11^^/ in the energy balance equation is the convolution of two rank- two 



jpp' • -^pp 
tensors and its expression in the thin-disk approximation (neglecting the off-diagonal components) is as follows 



(Vt^)w : n. 
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